Binary Inspirals in Nordstrom's Second Theory 
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We investigate Nordstrom's second theory of gravitation, with a focus on utilizing it as a testbed 
for developing techniques in numerical relativity. Numerical simulations of inspiraling compact star 
binaries are performed for this theory, and compared to the predictions of semi-analytic calculations 
(which are similar to Peters and Mathews' results for GR). The simulations are based on a co- 
rotating spherical coordinate system, where both finite difference and pseudo-spectral methods are 
used. We also adopt the "Hydro without Hydro" approximation, and the Weak Radiation Reaction 
approximation when the orbital motion is quasi-circular. We evolve a binary with quasi-circular 
initial data for hundreds of orbits and find that the resulting inspiral closely matches the 1/4 power 
law profile given by the semi-analytical calculations. We additionally find that an eccentric binary 
circularizes and precesses at the expected rates. The methods investigated thus provide a promising 
line of attack for the numerical modeling of long binary inspirals in general relativity. 



I. INTRODUCTION 

Numerical methods capable of modeling the strong 
gravitational fields of merging black holes and neutron 
stars were found several years ago after many years of 
searching. The main techniques are centered around the 
BSSN (see JlMd) and Generalized Harmonic (GH) for- 
mulations [y, [3 of general relativity. These two classes 
of techniques generally agree in their predictions 
and also match the Post-Newtonian calculations during 
the inspiral (101] . We are designing our own numerical rel- 
ativity techniques for comparable mass binaries with the 
goal of modeling long stretches of the inspiral phase pre- 
ceding the merger (as could be useful in discriminating 
neutron star equations of state [HI). We employ corotat- 
ing spherical coordinate systems to this end, as they are 
well adapted to the slow evolution of the binary during 
this phase. This is similar to the approach of the Caltech 
group ([H-[I3l) although we differ in most of the details. 

Spherical coordinate systems have several attractive 
features from a computational standpoint. A spherical 
mesh has a high density of grid points near the origin 
where it is needed to resolve the high curvature potential 
wells of the compact bodies. This thus provides an alter- 
native to the Adaptive Mesh Refinement methods (see 
e.g. [Ill, [III ) which have seen success lately. Addition- 
ally, a spherical grid provides a S*^ outer boundary, which 
is more suited to outgoing waves than a cubic grid. One 
of the key features is that it allows for the use of pseudo- 
spectral methods; by using spherical harmonics the 3-f 1 
problem can be split into a set of 1-1-1 problems which can 
be solved with fast and accurate implicit methods. The 
usage of spherical harmonics likewise avoids the compli- 
cations due to the coordinate singularities in 9 and (j)- A 
spherical system also meshes nicely with corotating co- 
ordinates, as proposed in [l^j- This removes most of the 
dynamical motion in the grid, thus helping to cut down 



on spurious field excitations. 

We are currently supplementing this primary numeri- 
cal method with some useful additional approximations. 
Neutron star like compact bodies are our focus, and thus 
the complications of modeling black holes are avoided. 
Furthermore we make use of the "Hydro- without-Hydro" 
approximation so that a relativistic fluid dynamics 
code is not currently needed. Radial profiles for isolated 
polytropic stars are generated and then inserted into a 
binary system. The gravitational fields and net accelera- 
tions of the stars in the binary are evolved, while the ra- 
dial density profiles of the stars are held constant. This is 
a good approximation as geodesic deviations due to tidal 
distortions of the stars only arise at high post-Newtonian 
order. 

Another useful approximation is possible when the 
stars are in quasi-circular orbit (which is generally the 
case for astrophysical binaries as the gravitational ra- 
diation circularizes the orbit). As pointed out in jl7| . 
in a corotating coordinate system a quasi-circular bi- 
nary evolves at the slow radiation reaction time scale 
T ^ d^/M^ for a separation d and mass M. This 
makes the Weak Radiation Reaction (WRR) approxima- 
tion possible; the second time derivatives in the field 
equations can be dropped, thus changing the character 
of the differential equations. 

There are some numerical drawbacks to using spherical 
coordinate systems. A primary one is that spherical har- 
monic decomposition and synthesis are computationally 
expensive. For instance the decomposition; 



Ai,n{r) = / f{r,9,^)Yt:,{d,cj))dn 
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is an 0{N^) calculation if discretized in a direct man- 
ner, although this can be reduced to 0{N'^) by defining 
intermediate variables and splitting the integration into 
two steps. The computational cost is further reduced in 
our case as we only need decompose and synthesize in 
the small volumes containing our compact bodies. This 
makes it possible to run high resolution meshes in a rea- 
sonable amount of time on a single computer. Further 
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reduction the the order of these operations is also possi- 
ble: [l9l - l2]| , although these methods have a large amount 
of overhead and only become efficient for large values of 
I and m. 

There are some other subtle issues that arise in spher- 
ical coordinate systems. For instance when using second 
order accurate methods one finds that very fine meshes 
are needed to accurately compute the accelerations of the 
stars. We find that using higher order methods allows for 
accurate results with much lower resolution. This is de- 
scribed in greater detail in section 3. 

Before we attempt to use these techniques in a full 
GR code we have elected to test them in a scalar gravity 
model theory (following the advice of |22j] and |2S]). This 
forms the focus of this paper. After the tests are success- 
ful we can then feel confident in applying the methods to 
a full general relativistic code. After considering several 
scalar theories, we chose Nordstrom's second theory (see 
[23|V It is a fully conservative metric theory, with param- 
eterized post- Newtonian parameters 7 = — l,/3 = l/2 
and all others zero. The theory thus has no Nordtvedt 
effect (see [25||), and the stars move on Keplerian orbits 
in the limit of small mutual gravitational potential. This 
is nontrivial since the Nordtvedt effect can produce con- 
siderable deviation from Keplerian orbits for highly com- 
pact bodies, even at arbitrarily large separations. We 
can therefore perform a semi-analytic calculation of the 
rate of energy and angular momentum loss for a binary 
based on the Keplerian motion, in a similar fashion to 
Peters and Mathews' calculation in GR (t2|j,[23). We 
find that the radiation is quadrupole to the leading order, 
despite being based on scalar field, due to approximate 
conservation of mass and momentum. 

The fully conservative nature of Nordstrom's theory 
has an additional benefit in that it justifies the use of the 
"Hydro-without-Hydro" approximation. The approxi- 
mation is safe since corrections to the star's equilibrium 
do not show up through first post-Newtonian order for 
fully conservative theories [28j . 

In the following sections we first review the details of 
Nordstrom's second theory. We describe the construction 
of isolated stars in this theory, and the equations of mo- 
tion for the matter. We then review the IPN equations 
of motion, which lead to precession for eccentric bina- 
ries. We calculate the energy and angular momentum 
loss on a large S'^ in the wave zone, and apply the results 
to a Keplerian system, which cause the binary to decay 
and circularize in a similar fashion to GR. Next we de- 
tail the numerical methods used to simulate the binary, 
including switching to a co-rotating reference frame, and 
several preliminary test cases that are needed to check 
the performance of the code. We finish with the results 
of simulating Nordstrom binaries in both quasi-circular 
and slightly eccentric orbits. 

We obtain a nice agreement between our simulation 
and the analytic calculations. The orbital decays of 
quasi-circular and eccentric binaries closely match the 
expectations from our semi-analytic calculations, as seen 



in the evolution of the semi-major axis a and eccentricity 
e. We also find that the eccentric binaries precess at the 
rate predicted by 1st Post Newtonian calculations, which 
gives supporting evidence that Nordstrom's theory obeys 
the Strong Equivalence Principle (SEP). 



II. NORDSTROM'S SECOND THEORY 

After the development of special relativity it became 
clear that Newton's gravitational theory could no longer 
be completely correct. For instance the Poisson equation 
is solved simultaneously throughout all of space, and thus 
the acceleration of a massive body would allow for the in- 
stantaneous transmission of signals, in disagreement with 
the finite speed of light. Comparisons with electromag- 
netic theory suggested that the spatial Laplacian oper- 
ator should be replaced with the D'Alembertian wave 
operator. Gunnar Nordstrom used this idea to develop 
two relativistic scalar gravity theories a couple years be- 
fore Einstein discovered General Relativity (we focus on 
the second, more complete theory). Nordstrom's the- 
ory is not a viable candidate for describing relativistic 
gravitation, as Nordstrom quickly realized after Einstein 
presented his tensor theory. For instance, being confor- 
mally flat, it predicts zero bending of light. However, it 
does provide an excellent arena for developing tools for 
numerical relativity. 

Einstein and Fokker demonstrated that Nordstrom's 
second theory could be expressed in geometric terms. 
The metric is conformally flat (thus the Weyl tensor is 
zero: Cap-yS = 0) and is specified by a scalar field tp: 



9iJ.y = (1 + vfVt^i' 



(2) 



The scalar field is generated in turn by setting the Ricci 
scalar to be equal to the trace of the stress energy tensor 
(see [H]): 



R ^ 247rT. 



(3) 



For our simulations we will use the standard perfect fluid 
stress energy tensor: 

T^"^ = {p + pe + p)u^u'' + g'^y, (4) 

so the held equation becomes: 

□(^ = 47r(l + (^)3(p + pe-3p). (5) 

Being conformally flat, Nordstrom's theory effectively 
has a background geometry rj^i, , and a conserved energy- 
momentum complex can be constructed: 



I 

47r 



+ (1 + </j)6t'^-. 



(6) 



with t^^'^ — 0. Nordstrom's theory has gravitational 
waves, although they are different from those in gen- 
eral relativity since Cap-yS = 0. Namely, the waves 
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occur in the scalar field Lp as described by equation 
([5]), and ^ describes their localizable energy density: 
tO" = (l/87r)[(9,^)2 + (V^)2]. 



A. Single Star Solution 

We utilize the hydro-without-hydro approximation in 
our code. This is a good approximation since tidal distor- 
tions of the stars do not affect the orbital dynamics until 
relatively high Post-Newtonian order. We thus need to 
construct a model of a single star in Nordstrom's theory 
to use as a source in the main code. For a single isolated 
star equation ^ becomes: 

^dr{r''dr^) = 47r(l + ^f[p{l + e) - 3p]. (7) 

We will assume the star is a polytrope, so that the pres- 
sure and internal energy are given by: 



r - r 



(8) 
(9) 



which reduces equation ([7]) to two variables. Applica- 
tion of the projection operator = u°'u^ + 5'^ on the 
divergence of the stress energy tensor T'^'^-f^ = gives: 



(10) 



where h is the relativistic specific enthalpy: h = h{p) — 
l + e+p/ p. With w° = 1/(1 + (y5) and = this reduces 
this to: 



dtp 
dr 



(l + (^) dh 
h dr 



(11) 



We solve ([7]) and (jTT]) iteratively. In the first iteration 
the (1 + ipY Slid (1 + If) terms are dropped respectively, 
resulting in a modified Lane-Emden equation which is 
solved by specifying the central pressure pc and integrat- 
ing out. For successive iterations the {\ + ipY and i^ + f) 
terms are reinserted with the previous iteration's solution 
for if. The iterative process converges for values of pc be- 
low some critical value (which depends on k, and F). For 
instance, with a stiff equation of state with F = 2 we 
find the most compact star that can be formed has ra- 
dius R ~ 1.75M, which is smaller than the event horizon 
radius in classical GR (in addition to being smaller than 
the Buchdahl-Bondi bound of i? = 9/4M for GR). We 
use the softer value of F = 5/3 to generate stars for the 
simulation, as the smoother transition to zero density at 
the outer boundary of the star allows the spherical har- 
monic decomposition to converge more quickly. 



B. Matter Equations of Motion 

The isolated star solutions will be used in the binary 
simulation with the hydro- without-hydro approximation. 



so that only the net accelerations of the centers of masses 
need to be found. The conserved density p* (see |2§|) 
proves useful to this end: 



I ^l/2 i'^ + vf 
P = P[-9) ' u = p 

vi — 



(12) 



with dtp* + di(p*v'^) = 0. The conserved density can be 
used to define the center of mass for star a: 



1 



p*x^d^x 



(13) 



where M* is volume integral of p*. Differentiating twice 
gives the coordinate acceleration ajj of the center of mass: 

ai^^ f p*^d^x. (14) 

An expression for the local accelerations dv^ /dt is then 
found by expanding the projection of the divergence of 
the stress energy tensor ([TU)) : 

-(«°)-24q-'>,.. (15) 
pn 

During the simulation the terms in (1151) are calculated 
and integrated with (IT4l) to find the net acceleration of 
the star (with the radial profile for p* kept constant). 



Semi- Analytical Calculation of Orbital 
Evolution 



Semi-analytic calculations of the evolution of a binary 
system in Nordstrom's theory are needed so that com- 
parisons can be made with the numerical simulations. 
We first review the IPN results for the theory, which 
follow from having the parameterized post-Newtonian 
parameters given by 7 = —1, /? = 1/2 (and all others 
zero). The theory therefore has no Nordtvedt effect (since 
4/3 — 7 — 3 = 0- see [2^) and compact stars move on 
Keplerian orbits at large separations. The IPN analysis 
also leads to precession for eccentric binaries. Next we 
will solve for the energy and angular momentum radiated 
at large radius by a binary on a Keplerian orbit, and then 
apply these results to evolve the orbital parameters. 

Consider first a binary with stars a and 6 in a slightly 
eccentric orbit with an angular velocity w. Arrange a 
Cartesian coordinate system so that the binary is instan- 
taneously aligned along the x-axis, with a separation d 
and the bulk of the velocity in the y direction: ^ ujd/2, 
and a small amount of radial velocity v^. The PN anal- 
ysis then finds that the acceleration Xa of star a in the x 
direction is: 



Mb 

' d2 



1 



Mb 
d 



<-7^{vl?-{vl-vlH 



(16) 
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where Mb = J^p*{l + {l/2)v'^ - {1/2)U + e)d^x is the 
gravitational mass (equal to the inertial mass) of star b 
(with U — /jP'|x — x'\~^cfix' and v = v — vq where vq 
is velocity of the center of mass). The acceleration in the 
y direction is: 



Va = 



Mbj yy - vl)vt 



(17) 



At large separations the post-Newtonian corrections 
Mjj/d and drop out, leaving the binary in a Keplerian 
orbit. This confirms that Nordstrom's second theory has 
no Nordtvedt effect through IPN order (and as a purely 
geometric theory it should follow the SEP). One further 
consequence of the post- Newtonian corrections ([T5)) and 
(|17p is that an eccentric binary will precess. A binary 
with semi-major axis a, eccentricity e, and total mass M 
will precess Aw radians per orbit [29]: 



Aw 



a(l-e2)' 



(18) 



This is six times slower and in the opposite direction as 
GR. 



D. Energy Loss at Outer Boundary 



This can be solved with a Green's function: 

[p*'(l - (l/2)«'2 + e'- ip'/p' + 0{v^))Ut ^3 , 



J |x-x'| 

(25) 

which is evaluated at the retarded time t' = t — |x — x'|. 
The l/|x — x'l term expands into: 



1 



1 x^ x^ 



+ 



(26) 



Only the first term is needed since it is utilized at large 
radius. Likewise only the first term in the expansion: 



r — X — X = 



L^x^' x^x^ {x^' x^' — 



(27) 



will be needed. Equation (^5]) can thereby be expanded 
out in multipole moments: 



if : 



(28) 



with 

M = 



J p*'il~{l/2)v'^ + e'-3p'/p' + 0{v*))d^x', (29) 



We next find the energy radiated by a binary due to the 
waves generated in ip. The energy-momentum complex[6] 
follows the standard fiat-space conservation law: 

t"-'^^ = 0. (19) 

Integration inside a volume bounded by S'^ gives: 

t°°,od'x^~ Jt'>\,d'x, (20) 

or, using Gauss's law: 

dtE^- j t°'n,dS. (21) 

with i*'* given by: 



t' 



Oi 



dt(pdi(p. 



(22) 



The spatial derivative diip can be transformed into 
—rfdtip since at large radius Lp is approximately a spher- 
ical wave: <p ~ sin(t — r)/r. This leads to: 



dtE = -^J {dt^fr^dn 



(23) 



for the energy loss. 

An expression for dt(p is needed next. Rewriting the 
equation of motion for the field ([5]) with the conserved 
mass density p* gives: 

□V5 = 47rp*(l-i;2)i/2(l + e-3p/p). (24) 



D' = I p*\l - (l/2)w'2 +e'- 3p'/p' + 0{v^))x''d^-' 



(30) 



and 



Q«J= / p*'{l-{l/2)v''- + e' -3p'/p' + 0{v'^))x"x^'d^x' 

(31) 

(and higher multipoles unnecessary for a leading order 
calculation of the energy loss). 

The time derivatives of these multipole moments are 
needed for 



The quadrupole contribution SfQ* 
scales as v^, so all terms of order and higher will 
be dropped for this leading order calculation. We first 
examine the time derivative of the monopole: 



dtM- 



p*'-{l - {l/2)v'^ + e' - ?,p'/p')d^x', (32) 



where {d/dt) Jy p* fd^x = Jy p*{df /dt)d^x has been used 
to pass the time derivative through p*. Note that if 
the two stars are on a quasi-circular orbit then the time 
derivatives in the integrand are on the radiation reaction 
timescale, and thus contribute radiation at a far lower 
scale than v^. However, they do contribute to eccen- 
tric binaries where the separation oscillates on the orbital 
timescale. 

First consider the de/dt term in ([5^ . and to sim- 
plify assume that F = 2, although the result holds in 
general. The specific internal energy can be expanded: 
e = Kp = Kp*{l - (1/2)^2 - 3(p + 0{v^)), where only the 
first term can contribute through quadrupole order. This 
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term reduces to de/dt = —Kp*diV^. Thus any monopole 
radiation of order from the de/dt term stems from the 
"breathing" motion as the star expands and contracts 
during the elhptical orbit. However, Nordstrom's sec- 
ond theory is fully conservative, and accordingly the stel- 
lar matter undergoes no "breathing" motion: the central 
density is constant to IPN order (which justifies the use 
of the "Hydro- without-Hydro" assumption). Therefore 
the entire de/dt term does not contribute at order. 
The same holds for the d{3p/p)/dt term. We find the 
monopole contribution to the radiation; 



this becomes: 

(39) 

This reduces to: 



(40) 



for zero eccentricity case, which is six times smaller than 
the value found in GR. 



dtM 



p*'dt{v'^)d^x' + 0{v^), (33) 



which is at order if the orbit is eccentric and is much 
smaller otherwise (note also that the total time deriva- 
tive has been switched to a partial derivative since the 
velocity is now essentially constant throughout the star). 
Note that it is also convenient that radial pulsations do 
not contribute at leading order as this allows the stars to 
be treated as point bodies in later calculations. 

We next examine the contribution from the dipole via 
dfD^. Distributing the two time derivatives through the 
integrand gives: 

d^D^ = I p*'[{l - v'^ +e'- Sp'/P')^ (34) 

dt \ dt J dt ' 

The / p*' {dv^' /dt)d^x' term integrates to zero due to 
Newton's third law. Any corrections would enter at 
order at the soonest, and all the other terms in equation 
(p4|l also scale as . Thus the dipole does not radiate to 
leading order: 



(35) 



The final multipole moment that can contribute 
through leading order is the quadrupole. The (l/2)w'^, 
e', and — 3p'/p' terms all enter at 0{v'), leaving only: 



dfQ'i = / p*'x"x^'d^x' + 0{v'^). 



(36) 



Combining the terms from the monopole and 
quadrupole gives: 

dtE^~^ J {dtM + {l/2)n,n,dfQ'^fdn (37) 

for the rate of energy loss. With the integrals: 

f Att ■ 

ninjdn = Y^''^ (38) 



E. Angular Momentum Loss at Outer Boundary 

A similar calculation determines the rate at which the 
angular momentum is radiated. With equation (|19l) and: 



U^e,jk J xH°''d^x 



(41) 



fe / xH^\id^x = -eijk I xH^^nUS. (42) 



one finds: 

dU ^ 
dt 



The spatial stress energy tensor components are: 

t'^' = ^dk^dw - ^<5f r/"^a„^5^^, (43) 

47r OTT 

and thus: 

^ = J x^dkfdiipn'dS. (44) 

The ri^^daipdpip term is zero due to the anti-symmetry 
of eijk- 

Only one spatial partial derivative should both con- 
verted into a time derivative via diip = —n^dtif- 



^ = -^e^/c j x^dk^^idtM + {l/2)ninmd^^Q'^)dS, 

(45) 

as converting both yields zero due to the anti-symmetry 
of eijk (note also that the nid^D'- term has been dropped 
from dtf since it doesn't contribute at leading order). In 
particular, examination of eijk shows that the dkf term 
that needs to be expanded to higher order, while the ditp 
can approximated by the time derivative. Application of 
the spatial derivative dk to the monopole term in ([28)1 
gives: 



5^ ( ~-/n - 



(46) 



Insertion of this term into (|45p also gives zero due to 
the anti-symmetry of tijk, and the same applies to all 
terms stemming from derivatives of powers of r. The 
only remaining relevant terms in dk^ are: 



(47) 
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Upon expansion in (1451) the —dtD^/r^ terms are multi- 
plied by either one or three v} terms, and as J nidU, = 
J ninjUkdn — these dipole terms integrate to zero. 
Equation (|45| thus becomes: 



(48) 

Reusing the integrals ([551) leads to: 



(49) 



where the dtM term has also dropped out, again due to 
djk- This expression is precisely one sixth of the value 
given by general relativity. 



F. Application to Keplerian Orbits 

We can now apply the equations for the rates of energy 
(p9)) and angular momentum loss p9)) to a binary star 
system in Keplerian orbit. The stars have gravitational 
masses mi and TO2, a semi-major axis a and eccentricity 
e. The separation d between the two stars is determined 
by the phase 0: 



a(l-e^) 
1 + ecos((/)) ' 



(50) 



with the distances of the stars from the center of mass 
being: 



di 



m2 



mi + 7712 



d, 



mi 



mi + 7712 



d. (51) 



The non-zero quadrupole moment components are: 

Qxx = fJ-d^ cos^ 0, 
Qyy — ^d^ sin^ (j), 
Qxy = Qyx = ^^d^ sin 4> cos (j) 



(52) 



with /i — 77117772/(7771 + m2). Note that the Qij are ex- 
pressed here in terms of the gravitational mass instead of 
the conserved mass 777* = J p*d^x as used in the previous 
section. Switching between the two involves corrections 
of order which do not affect the leading order calcula- 
tion of the orbital evolution. The second and third time 
derivatives of the quadrupoles are needed for (p9| and 
Using the angular velocity: 



[{mi + 7772)a(l - e2)]i/2 

d2 



(53) 



one finds the second derivatives to be: 
d Qxx 



dt^ 



7(4cos(20) -I- e(3cos(0) -f cos(30))), 



vv 



dt^ 



= 7(4cos(20) + e(4e + 7cos((/)) + cos(3(/)))), 



^ Qxy d Qyx 

~df~ " ~dF~ " 

where 7 defined as: 



-27sin((?!))(4cos((/)) + e(3 + cos(20))), 

(54) 



77117772 



' 2a(l - 
The third derivatives evaluate to: 



(55) 



d^Q. 



dP 



— — ^ = I3{1 + e cos 0)^(2 sin 20 + 3esin0cos^ 0), 
dt-^ 

yy = -/3(1 + e cos (j)f{2 sin 20 + e sin 0(1 + 3 cos^ 0)), 



d Qxy d Qy 



dt^ dt^ 

— /3(1 -I- e cos 0)^(2 cos 20 — e cos 0(1 — 3 cos^ 0)), 

(56) 

with /3 defined as: 



2 Am\m\[mi+m2) 



(57) 



The time derivative of the monopole dtM is also needed 
for eccentric binaries: 

dtM — —^{midtVi + 77729(f|) = ^/3esin0(l -I- ecos0)^. 

(58) 

The energy therefore radiates at the rate: 

dtE ^ -— /32(l-|-ecos0)'*(4 + 2e2 + 8ecos0-|-2e2 cos^ 0), 
15 

(59) 

and likewise for the angular momentum: 
2 

dtL, = /37(l+e cos 0)3(8-26^+126 cos 0-|-6e2cos2 0). 

15 

(60) 

In general relativity the gravitational wave energy can't 
be localized at individual points in space, and therefore 
the expressions for the dE/dt and dLz/dt need to be 
averaged over an orbit. While the energy can be localized 
in Nordstrom's theory, we will also average (|59l) and ((60)) 
in order to directly compare with general relativity. This 
is not a bad approximation since the orbital parameters 
evolve even more slowly than in GR. We find: 



16 m\m\{mi + 7772) 

15 a5(l_g2)7/2 



1 



13 2^ 7 4 

6 H 6 

4 16 



(61) 



which is similar to the value given by Peters [27| for gen- 
eral relativity (differing only in the coefficients): 



{E) = - 



32 m\m\{mi + 7772) 

"5" a5(l- 62)7/2 



73 

24' 



37 
96^ 



(62) 
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Likewise for the angular momentum we average to find: 



(4) = 



16 m\m\{mi + 1112)^^'^ 



15 



a7/2(i_e2)2 



(63) 



which again is one sixth the value found in general rela- 
tivity. 

Finally, by using 

a = -mim2/2E, (64) 
= TO^TO2a(l — e^)/(mi + m2) (65) 

(E) and {L^} can be converted into rates of change for 
the semi-major axis (d) and eccentricity (e): 



/ da 
and 



32 mir7i2(TOi + TO2) 
15 a3(l-e2)7/2 

18 TOim2(TOl -|- TO2) 

~T a4(l-e2)5/2 



13 

T 



1 + ^6^ + 



7 

16' 



1+18^ 



(66) 



(67) 



In the zero eccentricity case the semi-major axis inte- 
grates to: 



128 

a{t) — ( a''(0) —mim2{mi + m2)t 



1/4 



(68) 



Reproducing this 1/4 power law inspiral, similar to the 
expression found by Peters and Mathews for GR, is a 
major test of the Nordstrom simulation. 

We also simulate binaries that retain some eccentricity 
e. Equations (1551) and (|57| can be combined to directly 
compare how a and e evolve: 



16 



1 -I- ilg^ + J-e 



7 „4 



27e(l-e2) 1 + 



(69) 



An eccentric Nordstrom binary will therefore circularize 
as it evolves. Reproducing this rate of circularization, 
and the rate of precession given in (|18|) . are the major 
tests of eccentric binaries in our numerical simulation. 



matches the binary's angular velocity. The differential 
equation becomes (dropping the bar notation): 

-QjV + '^d^ip + 2ujdtd^ip - ui'^djip + V'^ip = (70) 
AT:{l + ipfip{l+e)-3p). 

Spherical harmonics are used to split this 3-1-1 differ- 
ential equation into a set of 1-1-1 radial equations, one for 
each I and m term: 

-dt^Vhn + imujipi„i + 2iTnu>dt(pi„i + ■m'^uj'^ipi.m + 

^ -2 = 4^-5;™ (i,r), (71) 

where the source term is also decomposed: 
SiUt,r)= f dnYil,{9,cbm + ^f{p{l+e)~3p)]. (72) 



We have written a spherical harmonics package as part 
of the code in order to evaluate ([7^ and also synthesize 
the field (p from its (pim components. We will use ([7T|) to 
simulate binaries that retain some degree of eccentricity. 

For binaries that are in quasi-circular orbits, we make 
the Weak Radiation Reaction (WRR) approximation, 
which is based on the presence of two timescales: the fast 
orbital motion characterized by uj (with period T ~ i^~^), 
and the longer radiation reaction timescale t. For a bi- 
nary with separation d the ratio of these timescales goes 
as t/T - {d/M f/'^. Thus there is a hierarchy of scales 
among the time derivatives: 

d I d . I 

The second time derivative term is the smallest, and is 
dropped when making the WRR approximation, giving a 
first order in time differential equation (which resembles 
the 1-D Schrodinger equation): 

i dr[r'^dripim] ,,22,- • 

Ot^Piryi = ■7^z~\ 2 h (m w 4- imuj 



2rnuj 



-^l±^)^i„,^A7:SUt,r)] (74) 



III. NUMERICAL METHODS 

A. Co-rotating Coordinates and the Weak 
Radiation Reaction Approximation 

The field equation ([5]) needs to be finite differenced 
to simulate a binary system in Nordstrom's theory. We 
begin by rewriting the wave equation —dfip + S/'^ip in a 
corotating spherical coordinate system. The corotating 
azimuthal angle (j) is given in terms of the original by: 
= 0- n{t) (with n{t) = /q w(t')dt' and t, r, and 
remaining unchanged) where uj{t) always instantaneously 



for the TO 7^ terms and 

drir^dr^lo] 1(1 + 1) 



(^/0"4^5io(i,r-) =0 (75) 



for the time independent m — terms. 

To solve for the initial data for the Cauchy evolu- 
tion the single time derivative terms in equation (j74p are 
dropped, giving: 

dr[r^dripim] , / 2 2 + ^)\ 

+ [m CO -2 I (pim 

-47:Si^(t,r)^0. (76) 
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This leaves the boundary conditions. The inner bound- 
ary is given by setting the radial derivative of the 
terms to zero: 



0. 



(77) 



At the outer boundary the Sommerfeld outgoing wave 
boundary condition is used, which has the following form 
in the co-rotating reference frame: 



(78) 



The binaries are first modeled in a quasi-circular orbit 
(with the expected inspiral given by (j68p ) with equations 
([74| and (|75p used to evolve the field (p. As noted equa- 
tion (fM)) closely resembles the Schrodinger equation, so 
it is finite differenced using the Crank-Nicholson method 
(see e.g. 

The hyperbolic equation ((7T|) needs to be finite differ- 
enced for inspirals that retain some eccentricity. Given 
the success we have in evolving (TMl) with Crank- 
Nicholson, we devise a modified Crank-Nicholson for 
(ED: 



N+l 



"Pi 



2 ^ 



2^f -I- <^f + zm^At(<^f +1 - -^-1 



N+l 
i+1 



(79) 



Ar 



rAr 



,1(1 + I) 



2 2 
TO U! 



i+l 



l{l + 1) 2 2 ■ ■ \ N-l\ 

+ ( — ^ m LJ - tmojjip^ ) 

+2TTAt\S^+^ + Sf-^). 

We find that this balanced implicit method also allows 
for fast, stable, and accurate field evolutions. 



B. Numerical Tests 

We use a mixture of finite difference and pseudo- 
spectral methods, and these need to be tested first before 
evolving the full system. A useful test case is provided 
by the Newtonian system, with the inhomogeneous wave 
equation reducing to Poisson's equation: 



W LP = iirp, 
and the acceleration is given by = ma: 



(80) 



(81) 



Equation (|80p can be split into its spherical harmonic 
components by dropping the rn^uP'tpim term from equa- 
tion (Uni): 



^i™ + 47r5i™(i,r) =0, (82) 



1 1 — 

2.2 log(N) + 4.4 

"log_error.dat" 




FIG. 1: Plotted is the logio of the error e = sup |vJVetuton/v~ 
1| against the logio of the number of radial grid points A''. 
A linear fit gives a slope of —2.2 ± 0.2, i.e. second order 
convergence. The lowest resolution grid uses A'^ = 500 radial 
grid-points and L and M values up to 25. This is then doubled 
four times up to resolution with N = 8000 radial grid points 
and L and M values up to 400. 



with the source terms given by Sim{t,r) — 
f dnY*^{e,(j))p{t,r,e,(l}). Equation ^ is then fi- 
nite differenced in the standard second order manner. In 
order to check the convergence of the numerical solution 
we choose a matter distribution for the stars given by 
p{r) — pq{1 ~ / R^) (with some central density po and 
radius R) which allows for an exact analytical solution 
for ip. After solving for the (^/m components based on 
this source and then synthesizing the field we find 
that our numerical solution converges at second order to 
the analytical solution, as shown in figure ([T]). 

The Newtonian system presents a second test in ad- 
dition to the convergence of the field: the integral in 
equation ([ST|) must also converge to the expected in- 
verse square law. Note also that we are already using a 
"hydro-without-hydro" type of approximation: the mat- 
ter is described by a rigid density profile, so that only 
the acceleration of the center of mass needs to be de- 
termined. We find that accurately resolving this accel- 
eration is somewhat involved in a spherical coordinate 
system. For a given resolution, and using second order 
methods, we generally find that the solution for the field 
([50)1 is much more accurate than the result for the ac- 
celeration ((8T|) . The error in the computed acceleration 
also grows quite quickly as the separation increases. The 
error does converge to zero as the resolution is increased, 
but inconveniently high resolutions are needed to reduce 
it to acceptable values. 

This inconvenience arises because ([5T|) includes the 
gradient of the star's self field. The Newtonian self force 
integrates to zero analytically, but numerically there will 
be some small residue. Furthermore, the amplitude of 
the self field is much larger the I /r contribution from the 
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10 15 20 25 30 35 40 45 50 55 60 
Separation (d) 

FIG. 2: Plotted is the error = (a,.)/(l/d^) - 1 of the radial 
acceleration (compared to the expected Newtonian inverse 
square law with a mass of 1) as a function of the separa- 
tion d for the 2nd order and 12th order derivative methods. 
The standard 2nd order accurate method gives an error of 25 
percent for a separation of 60M (for a particular grid reso- 
lution), while the 12th order method gives an error of ~ 0.3 
percent, which is the same accuracy the field is resolved to at 
this separation. 



other star, so a small percent error in calculating its gra- 
dient can swamp the correct contribution from the 
other star, especially as the separation increases and the 
resolution inside the stars decreases (for a fixed spherical 
grid). 

Experimentation shows that the error grows almost as 
quickly if the analytic solution for (p is inserted into the 
mesh and then finite differenced to calculate (|5T|). How- 
ever, if the analytical derivatives ditp are inserted into 
([5T|) then the result for the acceleration is much more ac- 
curate - now comparable to the accuracy of (|80|) . We thus 
need much more accurate derivatives, which are available 
through higher order methods. Experimentally, each in- 
crease in the order of the derivatives gave better results 
(as compared to the analytical solutions), so we used 12th 
order accurate derivatives, given by: 



diV - 27720Ag; ^^'^''^ ~ '^'^'^'-'^ 495(/?,_4 
-2200(^i_3 + 7425(^i_2 - 23760(^i_i 
-h23760c/j,+i - 7A25ip,+2 + 2200(p,+3 
-495<y9i+4 + 72(^1+5 - 5ipi+e]. 



(83) 



The resulting improvement in the acceleration is shown in 
figure ([2]). Note that we do not claim to have a 12th order 
accurate code, as the field ip is still solved for with 2nd 
order accurate implicit methods. However, these higher 
order derivatives do allow the net accelerations computed 
in ((5T|) to have the same level of accuracy as the solution 
for if. 

The next test is to promote the Laplacian operator in 



vy_dot(d) = 2/(15 d^(9/2)) 
"vy_dot.dat" 
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FIG. 3: Plot of the measured acceleration dv^/dt (i.e. in the 
(j) direction) as a function of separation, which closely matches 
the expectation given in 1851 



(84) 



equation (|80)) to a wave operator: 
□(y9 = Anp. 



The spherical harmonic decomposition for this in coro- 
tating coordinates was given in (j7ip and the initial data is 
found by solving ([7S)) . The homogeneous version of ([75)) 
is solved by the Bessel and Neumann functions ji(mujr) 
and r]i(mujr), or in the case of the Sommerfeld bound- 
ary condition (1751) . by an outgoing wave spherical Han- 
kel function. At large r this asymptotically approaches 
a sinusoidal function of frequency muj with a l/r enve- 
lope. Our finite difference code correctly reproduces this 
asymptotic behavior. 

This simple wave equation system also provides an- 
other test of the net acceleration of the center of mass. 
For this test we retain the Newtonian equation ([5T|) to 
calculate accelerations based on the field ip. In New- 
tonian language, the wave equation Dip gives rise to a 
radiation reaction potential $(''<=''c*) (when driven by the 
orbital motion of the compact bodies). This potential 
^{react) j-gg^j^g jjj ^ drag forcc acting on the stars, caus- 
ing the orbit to decay in a manner that closely matches 
the full theory. A quadrupole calculation similar to the 
full Nordstrom theory version finds that the orbit decays 
at the rate: d = — (64/15)m^/£i^ where d is the separa- 
tion and nia — mi, is the mass of the equal mass stars. 
Let the two stars be instantaneously aligned along the 
X-axis in a Cartesian plane, at ±d/2 away from the ori- 
gin, with Newtonian velocities = ±y{ma/2dy^^ . The 
gradient of then gives an acceleration opposite 

to the velocity (i.e., in the (j) direction): 



yy = (l/15)mI/V(d/2)9/^ 



(85) 



thereby causing a quasi-circular inspiral. The code accu- 
rately reproduces the inverse nine-halves scaling of this 
acceleration, as shown in Figure ([3]), where we have set 
nia = 1/2 (so that the total mass is 1). 
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FIG. 4: Separation d as a function of time for a quasi-circular 
inspiral. The inspiral closely matches the 1/4 power law result 
given in ((68)) . 



FIG. 5: Plot of the rate of energy loss of the binary dE/dt 
as a function of the separation d, compared to the theoretical 
value given in (|61|l . where ml = m2 = 0.5 and e = 0. 



IV. MODELING NORDSTROM'S THEORY 

With the initial tests passed we can now simulate 
Nordstrom theory. Again we solve for the initial data 
with (j76p , but we do not decompose the source as we did 
in (j72p . Instead we rewrite the source term using the 
conserved density p*, finding: 



SUt.r)^ j dnY;„,{e, cj,)[Vl - v^P*{1 + e - 3p/p)]. 

(86) 

This removes much of the nonlinearity from the equation 
of motion for the field. We find e and Sp/p using p* 
(which is held constant under the hydro-without-hydro 
approximation), and the value of ip. 

We first model quasi-circular inspirals by making use of 
the WRR approximation. Initial data is found by start- 
ing with a circular Keplerian binary, and then iteratively 
tuning the orbital parameters to reduce the initial ec- 
centricity. The field is then evolved forward in time via 
([7^ , which is finite differenced with the Crank-Nicholson 
scheme. Note that it is crucial to also time average the 
source Sim{t,r) in the Crank- Nicholson scheme, as oth- 
erwise the source and field will be a half time step out of 
sync which results in spurious accelerations. 

After each time step the field ip{r, 9, 0) is synthesized 
from its spherical harmonic components and used to solve 
for the acceleration using a discretized combination of 
(ITi|) and p^ . As noted, 12th order finite differencing is 
used for the derivatives of ip so that the net acceleration 
is accurate. The components of equation (1151) are also ex- 
panded out in the corotating reference frame. The coro- 
tating coordinate transformation adds shift vector like 
terms to the metric which in turn leads to new terms in 
the sum of connection coefficients P-' ov"v^ . These ad- 

ap 

ditions effectively add a centripetal force, so that when 
the correct uj is found (via iterative testing) the stars will 
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Time (t) 

FIG. 6: Separation d as a function of time for a binary with 
an initial eccentricity of e 0.021. 



be nearly at rest in the co-rotating frame (expect for the 
slow merger due to the radiation reaction). 

The masses of the stars are set to ma = mt — 1/2 
for both the quasi-circular and eccentric simulations. An 
initial central density of p* ^ 3.72 x 10"'* and polytropic 
parameters F = 5/3 and k — 7.29 are picked to give a 
compactness of nia/Ra — 0.1. The stars are then rescaled 
in the simulation to give a total mass of 1, so that each 
star has a radius of 5M . The outer radial boundary is 
set to router = 560M, which places it in the wave zone 
for an initial separation of 31M. The radial grid spacing 
is Ar = 0.4M, and spherical harmonic modes up to / = 
m — 80 are used (in general the discretization is chosen 
so that there is comparable resolution in the angular and 
radial directions). As we use Crank- Nicholson to evolve 
(1741) for the quasi-circular case we can utilize a time step 
of At = 6.0M. This is much larger than the value given 
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FIG. 7: Plot of the Runge Lenz vector A as it evolves, which 
demonstrates that binary precesses about 2.5 times over the 
course of the inspiral. The circularization of the binary can 
also be seen as the magnitude of A drops by about a factor 
of two. 
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measured directly at the outer shell of the domain. We 
find a close match to the analytical prediction given in 
The energy contained in the volume of the com- 
putational domain is also tracked using ^ and found to 
decrease in sync with the energy flux through the outer 
boundary. The quasi-circular inspiral therefore closely 
matches the analytical expectations. 

Next we examine binaries that retain a small amount 
of eccentricity. The WRR approximation is dropped and 
the field is evolved with the modified Crank-Nicholson 
scheme given in ([7^ . We select a value of a; which gives 
the binary a small amount of eccentricity (e ~ 0.02) at 
an initial separation of a ^ 16A/. With this small value 
of e the binary decays only slightly faster than a quasi- 
circular binary with the same semi-major axis. The sys- 
tem circularizes noticeably over the evolution (finishing 
with e ^ 0.01) as can be seen in Fig. ([6]). 

We first compare the rate of precession to the analyti- 
cal profile given by Post Newtonian calculations. We 
find that the code matches this expectation closely: for 
the chosen initial data the binary precesses about Stt ra- 
dians over the course of dozens of orbits as the separation 
d decreases from 16M to WM. This precession can be 
seen by fitting the Runge-Lenz vector A to the binary's 
orbital parameters (and then averaging over an orbit to 
remove any orbital timescale artifacts). A points along 
the semi-major axis, and is proportional to the eccentric- 
ity. It therefore compactly displays both the precession 
and circularization of a binary, as shown in Figure ([7]). 

The gradual circularization of the binary can also be 
seen in Figure ([5]), where an extracted value for e(a) (for 
simulations with a range of time steps) is compared to 
the result found by integrating the analytical form in 
The slow secular evolution of e is quite sensitive to 
accumulating errors, and we find that the time step needs 
to be reduced to At ~ Ar/4 in order to get accurate 
results. 



FIG. 8: Plot of the eccentricity e as a function of the semi- 
major axis a for simulations with a time step At ranging from 
2.0 to 0.25 Ar vs. the analytical expectation from (|69|) . 

by the Courant-Friedrichs-Lewy (CFL) condition, but is 
sufficient to resolve the slow evolution of the system on 
the radiation reaction timescale. A large time step also 
proves to be stable for eccentric binaries, but in this case 
there is dynamical motion on the shorter orbital time 
scale, and experimentation shows that time steps similar 
to the CFL condition are needed in order to get accurate 
results. 

The binary separation as a function of time for the 
quasi-circular evolution is given in Figure Q ■ The initial 
separation is set to be 31M and the stars then inspiral 
to a separation of IIM (i.e. almost to merger) over the 
course of ~ 600 orbits. The inspiral closely matches the 
1/4 power law predicted by the semi-analytical cal- 
culation. In Figure ([5]) we plot the rate of energy loss as 



V. CONCLUSIONS 

The numerical techniques we have developed work well 
in simulating binary inspirals in Nordstrom's theory. The 
co-rotating spherical coordinate system is quite sensitive 
to the 2.5 PN order radiation reaction force, and is well 
suited to conserving angular momentum, allowing the for 
accurate evolution of inspirals over many orbits (^ 600 in 
the case of the quasi-circular inspiral) . The code demon- 
strates this, as it accurately matches the analytic pre- 
dictions we made for Nordstrom's theory, from the 1/4 
power law profile for a quasi-circular inspiral and the cor- 
responding rate of energy loss, to the rate of precession 
and circularization of an eccentric binary. 

As a side benefit we find that Nordstrom's second the- 
ory is a useful laboratory for developing computational 
techniques for numerical relativity. It is a fully relativis- 
tic metric theory, and yet has fairly simple field equations 
which can be put in a nearly linear form. Furthermore 
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the theory follows the SEP and thus stars move on Keple- 
rian orbits to lowest order, and there are no star-crushing 
effects. 

Our main goal is to use our corotating spherical coor- 
dinates framework to model the late inspirals of binaries 



in GR. We explore a modification of a generalized har- 
monic formulation of the Einstein equations to this end 
in 31]. Given the success we have had in applying it to 
Nordstrom binaries, we feel confident that it should allow 
for fast, long and accurate simulations in GR as well. 
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